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detector array is linear and outputs signals from the photodiodes 
mounted therein to a preprocessor for smoothing, correcting and filter- 
ing and subsequent processing to transform the signal from that pro- 
duced by an x-ray originating from a fan beam source, e.g., in a polar 
coordinate system, into the signal which would have been produced by a 
detector in an array on which a parallel beam is incident on a Cartesian 
coordinate system. The transformed data is converted to a gray scale va- 
lue for a picture element having a specific position in the Cartesian 
coordinate system and output to an appropriate display. Data is taken 
at each incremental angle as the beam source and detector array rotate 
around a target object The method of reconstructing this back projected 
image involves correcting and smoothing the output signals, scaling 
those corrected and smoother signals, and convolving the scaled signals 
into data characterizing the ray which is incident on each individual de- 
tector element into the equivalent intensity data had the incident ray ori- 
ginated from a parallel beam source. Also provided are methods of cor- 
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A 

METHOD AND APPARATUS FOR 
COMPUTING TOMOGRAPHIC SCANS 

BACKGROUND OF THE INVENTION 

Field of the Invention 

The present invention pertains generally to the use of 
a detachable x-ray detector array in combination with an 
existing radiation therapy x-ray simulator to produce 
computed tomographic (CT) image reconstructions. In more 
5 detail, the present invention relates to an apparatus for 

producing a CT scan from the width-collimated fan beam pro- 
duced by an existing x-ray simulator and to a method for 
transforming the data produced by the detector array of the 
simulator into a back projected image of the target object, 

10 More particularly, this invention pertains to a method for 

reconstructing x-ray projection data to form CT images from 
divergent x-ray beams. The method can utilize noisy x-ray 
projection data, for example, data resulting from poor x- 
ray power supplies, and correct for geometrical errors 

15 resulting from mechanical defects in the motion of the beam 

projector that cause the central projection ray to move 
about the center of the detector array as a function of the 
angle of beam projection. 
Introduction 

20 X-ray computed tomography (CT) is a technique for 

obtaining cross-sectional reconstructions of three dimen- 
sional objects using x-rays. In the simplest example of CT 
imaging, a narrow beam of penetrating x-rays is scanned 
across an object or patient in synchrony with a radiation 

25 detector on the opposite side of the patient. If the beam 

is monoenergetic or nearly so, the transmission of x-rays 
through the patient is given by the equation 

I - l 0 exp(-/ix) [1] 
where the patient is assumed to be a homogeneous medium 

30 with the attenuation coefficient \l. If the x-ray beam is 

intercepted by two regions with attenuation coefficients ^ x 
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and /x 2 and thicknesses x x and x 2 , the x-ray transmission is 
characterized as 

I = I 0 exp - (Mix x + n 2 x 2) E 2 3 

This formula is generalized to many (n) regions with 

5 different linear attenuation coefficients with the argument 

of the exponent 
n 

2 fl L X L = (1*!*! + 1*2*2 + • • • M n X n ) [3] 

i=l 

10 n 

I = SI 0 exp - (fx L x L ) [4] 
i=l 

Separate attenuation coefficients cannot be determined with 
a single transmission measurement because there are too 

15 many unknown values of p L in the equation- However, with 

multiple transmission measurements at different orienta- 
tions of the x-ray source and detector, the separate 
coefficients can be distinguished so that a cross-sectional 
display of coefficients is obtained across the plane of 

20 transmission measurements. By assigning gray levels (see 

below) to different ranges of attenuation coefficients, a 
display is obtained that represents various structures in 
the patient with different x-ray attenuation characteris- 
tics. This gray scale display of attenuation coefficients 

25 constitutes a CT image. 

The numbers computed by the reconstruction algorithm 
are not exact values of attenuation coefficients. Instead, 
they are integrals, termed CT numbers, which are related to 
attenuation coefficients. On most newer CT units, the CT 

30 numbers range from -1,000 for air to +1,000 for bone, with. 

the CT number for water set at 0. CT numbers normalized in 
this manner are termed Houns field units and provide a range 
of several CT numbers for a one percent (1%) change in 
attenuation coefficients. 

35 To portray the CT numbers as a gray scale visual 

display, a CRT is used. The CRT includes a contrast 
enhancement feature that superimposes the shades of gray 
available in the display device (i.e., the dynamic range of 
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the display) over the range of CT numbers of diagnostic 
interest. Control of image contrast with the contrast 
enhancement feature is essential in x-ray computed tomogra- 
phy because the electron density, and therefore the x-ray 
attenuation, is remarkably similar for most tissues of 
diagnostic interest. These electron densities vary from 
3.07 x 10 23 elec/cc for fat tissue to 5.59 x 10 23 elec/cc 
for the densest tissue; bone. Lung tissue has a much lower 
electron density, 0.83 x 10 23 elec/cc, because of the 
alveolar and bronchial spaces. 

The first CT systems were introduced in approximately 
1971 by the EMI Corporation of England. These early 
systems used an x-ray source mounted in a gantry with 
detectors. The patient was inserted between the x-ray 
source and the detectors and the joined x-ray source and 
detectors were moved about the patient to obtain projection 
rays through the patient. These values were fed to a 
computer which then reconstructed a cross sectional image 
of the plane through which the pencil beam of x-rays 
passed. During this translational scan of perhaps 40 cm in 
length, multiple (e.g., 160) measurements of the x-ray 
transmission were obtained. Next, the angular orientation 
of the scanning device was incremented one degree and a 
second translational scan of 160 transmission measurements 
was performed. This process was repeated at one degree 
increments through an arc of 180 degrees so that 28,800 x- 
ray transmission measurements were accumulated. Those 
measurements were then transmitted to a computer equipped 
with a mathematical algorithms for reconstructing an image 
of attenuation coefficients across the anatomical plane 
defined by the scanning x-ray beam. 

Although this approach yielded satisfactory images of 
stationary objects, considerable time (4-5 minutes) was 
required for data accumulation and the images were subject 
to motion blurring. Soon after the introduction of pencil 
beam scanners, fan-shaped x-ray beams were introduced so 
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that multiple measurements of x-ray transmission could be 
made simultaneously. Fan beam geometries, with increments 
of a few degrees for the different orientations (e,g., a 
30-degree fan beam and 10-degree angular increments) , 
reduced the scan time to 20-60 seconds and improved the 
image quality by reducing the effects of motion- Computed 
tomographic scanners with x-ray fan beam geometries and 
multiple radiation detectors constituted the second genera- 
tion of CT scanners. 

In late 1975, the third generation of CT scanner was 
introduced. These scanners eliminated the translational 
motion of previous scanners, using rotational motion of the 
x-ray tube and detector array or rotational motion of the 
x-ray tube within a stationary circular array of 600 or 
more detectors. With these scanners, data accumulation 
times as fast as two seconds are achievable. 

Both stationary and rotating anode x-ray tubes are 
used in CT scanners. Many of the translation-rotation CT 
scanners have an oil-cooled, stationary anode x-ray tube 
with a focal spot on the order of 2 X 16 mm. The limited 
output of these x-ray tubes necessitates a sampling time of 
about 5 msec for each measurement of x-ray transmission. 
This sampling time, together with the time required to move 
and rotate the source and detector, limits the speed with 
which data can be accumulated with CT units using translat- 
ional and rotational motion. 

To reduce the sampling time of 2-3 msec/ most fast- 
scan CT units use rotating-anode x-ray tubes, often with a? 
pulsed x-ray beam, to achieve higher x-ray outputs- Even 
with rotating-anode tubes, the heat-storage capacity of the 
anode may be exceeded if cooling periods are not observed 
between sets of successive images. 

After transmission through the patient, the x-ray beam 
is collimated to confine the transmission to a slice with 
a thickness of a few millimeters and to reduce scattered 
radiation to less than one percent (1%) of the primary beam 
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intensity. The height of the collimator defines the 
thickness of the CT slice. This height, when combined with 
the area of a single picture element (pixel) in the dis- 
play, defines the three-dimensional volume element (voxel) 
5 in the patient corresponding to the two-dimensional pixel 

of the display. A voxel encompassing a boundary between 
two tissue structures (e.g., muscle and bone) yields an 
attenuation coefficient of the two structures. This 
"partial volume artifact" may be reduced by narrowing the 

10 collimator to yield thinner slices. However, this approach 

reduces the intensity of the x-rays incident upon the 
detector and the detector signals are subject to greater 
statistical fluctuations, thus introducing more noise in 
the displayed image. 

15 To reduce the detector response time, all detectors 

used in CT scanning are operated in current rather than 
pulse mode. Also, rejection of scattered radiation is 
assigned to the detector collimator rather than to pulse 
height analyzers. Detectors for CT scanning are chosen on 

20 the basis of detection efficiency (greater than 50%) , short 

response time and stability of operation, and are either 
gas-filled ionization chambers or solid scintillation 
detectors. With any detector, the stability of response 
from one transmission measurement to the next is essential 

25 for the production of artifact-free reconstruction images. 

With a pure rotational source and detector geometry, for 
example, detector instability gives rise to ring-shaped 
artifacts in the image. Minimum energy dependence of the 
detector over the energy range for the CT x-ray beam also 

30 is important if corrections for beam hardening are to be 

applicable to all patient sizes and configurations. 

All of the early CT systems were designed and built 
only to perform CT studies. The concept of using other 
types of radiation sources that had not been specifically 

35 designed for CT imaging was initiated in the mid 1970 's. 
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Several of these efforts utilized existing x-ray- 
therapy simulators- An x-ray simulator is a device that 
duplicates a radiation treatment unit in terms of its 
geometric, mechanical and optical properties, but uses a 
diagnostic x-ray tube as the source of radiation to simu- 
late the properties of the treatment beam* A simulator 
allows the beam direction and the treatment fields to be 
determined while encompassing the target object with the 
simulator's irradiation. Since the simulator's emissions 
are generally less intense and less energetic than the 
emissions of therapy devices, there is a reduction in the 
target object's exposure to radiation. 

The combination of a detector system and an x-ray 
therapy simulator provides the necessary front end of a CT 
system. Application of the requisite information process- 
ing techniques and algorithmic reconstruction processes in 
combination with the simulator cum detector system enable 
production of CT images. 

In another type of radiography, known as panoramic 
radiography, the x-ray source and detector lie in the same 
plane just as is the case for an x-ray simulator, but they 
are rotated in a plane that is transverse to the plane 
containing the source and detector. In panoramic x-ray 
devices such as those used to produce images of the maxi- 
llo-facial region, the center of rotation of the source and 
detector moves along an elliptical path, the output signal 
of the detector being read periodically and then condi- 
tioned to produce a panoramic x-ray image of the target* 
object. In contrast with a CT scanner, such devices 
typically scan only approximately 230-240° about the target 
object. 

Reconstruction Algorithm 

The first reconstruction formula was derived by an 
Austrian mathematician, P.J. Radon. Uber die Bestimmung 
von Funktionen durch Ihre Integralwerte Laengs Gewisser 
Mannigfaltigkeiten, Berichte Sachsische Akadamie der 
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Wissenschaf fen (Leipzig) , Mathematische-Physische Klasse 
69 , 262-277 (1917) . His methods were based on the assump- 
tion of use of a parallel beam and concomitant data collec- 
tion. The fundamental process can be described as a 
5 convolution followed by a back projection. Subsequent 

methods derived from the Radon method utilized variations 
in the application of filters and convolution to correct 
for the effect of noise in the reconstruction process. 
Pavkovich and Nunan (U.S. Patents 4,149,248 and 

10 4,365,339) modified Radon's algorithm to accommodate fan 

beam geometry. Their method consisted of a convolution of 
the original projection data with a function derived from 
the inverse Fourier transform (IFT) of a product of the 
Fourier transform of the projection data and the ramp 

15 function in frequency space. These results were then back 

projected to the point of interest, (X,Y) . Their convolu- 
tion was essentially an angular convolution for any given 
point. All of their mathematical manipulations were 
performed in configuration (real) space with no resort to 

20 transforms to frequency space. 

Algorithms based upon the Radon method such as that 
developed by Pavkovich and Nunan suffer from two defects: 
(1) they are computationally burdensome, requiring both 
excessive memory and time of execution, and (2) they are 

25 very sensitive to noisy data; that is, resolution of the 

reconstructed image decreases disproportionately with 
increases in noise in the data. The Radon reconstruction 
method is theoretically perfect when . operating on perfect 
data; that is, an infinite number of parallel, noiseless 

30 projections detected by an infinite number of detectors of 

infinitesimal dimension. These requirements unfortunately 
do not exist in physical systems. Using real data, results 
from the Radon method for calculating attenuation coeffi- 
cients degenerate rapidly as the data becomes noisy and a 

35 finite number of discretely sized detectors is used. 
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Known reconstruction algorithms are of four general 
types: (1) simple back projection (Radon); (2) integral 
equations; (3) Fourier transform; and (4) series expan- 
sion. In the present invention, the Fourier transform is 
5 utilized. In this approach, the x-ray attenuation pattern 

at each angular orientation is separated into frequency 
components of various amplitudes, similar to the way a 
musical note can be divided into relative contributions of 
different frequencies. From these frequency components , 

10 the entire image is assembled in "frequency space" and then 

reconstructed by an inverse Fourier transform reconstruc- 
tion process into a spatially correct image. 

While the various algorithms used for CT image recon- 
struction each have their own limitations, the quality of 

15 the overall procedure is dominated by the quantity and 

quality of the measured transmission data. The quantity of 
data is restricted by the specific scanner design and by 
limitations placed on time and computer resources. Phenom- 
ena which tend to degrade the quality of the measured data 

20 include: 

a) Geometrical errors such as misalignment or motion 
of the scanning system or patient motion, 

b) Instability of the x-ray source, 

c) Statistical fluctuation of the measured signal, 
25 d) Polychromaticity (non-monochromaticity) of the x- 

ray beam, 

e) The finite dimensions of the scanning aperture, 

f ) Residual signal due to the time response function 
of the detector system (afterglow) . 

30 If, as a result of these factors, the projection values 

derived from the measured data do not adequately represent 
the line integrals of the linear attenuation coefficients 
within the slice being scanned, even the most perfect 
reconstruction algorithm will give rise to a distorted 

35 image. Each of these factors and the manner in which they 



WO 93/00649 PCT/US92/05276 

7 

are addressed in the case of the method of the present 
invention is set out in the following paragraphs. 
Geometrical Error 

The most practical approach to error introduced by 
5 scanner misalignment or motion lies in the construction and 

maintenance of the scanner itself, including periodic 
testing of the mechanical registration.. Patient motion is 
less controllable but .can usually be minimized by proper 
patient support and through the use of the fast (fan beam) 

10 scanner. If patient motion is monitored, such as through 

the use of a transducer arrangement or a laser beam reflec- 
tion method, then data correction is feasible. Using the 
method of the present invention, correction can even be 
made for eccentricity in the motion of the center of 

15 rotation; that is, when the center of rotation does not 

coincide with the center of reconstruction at a given 
range. Even such gross eccentricity as results from the 
use of a panoramic x-ray device is measured and subse- 
quently corrected for all angles involved in the back 

20 projection and reconstruction. 

X-Ray Source Instability 

The instability of the x-ray source is generally 
corrected for by adjusting the measured data in accordance 
with the signal measured by a reference detector. Another 

25 approach is to monitor the electronic parameters of the x- 

ray source, such as kVp and mA, and to make the corrections 
to the measured data from this information. The reference 
detector method, while not extremely sensitive to kilovolt- 
age variations, is simpler and more easily utilized through 

30 electronic hardware. Computer correction in accordance 

with the present invention makes possible the use of lower 
cost x-ray power supplies and kVp monitoring to correct for 
changes in effective beam energy. 
Statistical Fluctuation 

35 Systemic error may arise in such forms as drift and 

gain variations in the detector system and associated 
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electronics, or in the form of background or bias currents. 
If these variations can be monitored and quantified, 
corrections can be implemented through hardware circuitry 
or computer correction. 
5 Polychromatic Effects 

Since flux-rate requirements based on statistical 
considerations for any reasonable scan time normally rule 
out the use of other sources of radiation, an x-ray tube is 
the only practical photon source for CT scanning. As a 

10 result, a spectral distribution of photons is involved 

rather than photons of a single energy. Since lower energy 
photons have higher attenuation coefficients than higher 
energy ones, the beam becomes progressively "harder" as it 
traverses an increasing patient thickness. A "harder" 

15 beam, having a lower effective attenuation coefficient than 

a "softer" one, introduces a degree of inconsistency into 
the measured data used for reconstruction. In the absence 
of compensation for this effect, e.g., utilizing a fixed 
length water bath or software corrections applied at the 

20 preprocessing stage, this effect leads to a distorted image 

characterized by a general increase in coefficients from 
the center to the periphery of the cross-section. 

The object of preprocessing the measured transmission 
data before reconstruction is to "linearize" the logarithm 

25 of the ratio of the incident to transmitted intensity. For 

a monoenergetic beam traversing a homogeneous medium, this 
logarithm is a linear function of increasing thickness, 
while for a polychromatic beam this function is no longer 
linear. If the characteristic attenuation of a particular 

30 incident x-ray beam by an increasing mass thickness of 

water is known and if it can safely be assumed that the 
materials encountered within the body have attenuation 
properties similar to those of water, the measured data can 
be corrected to an idealized monoenergetic (linear) re- 

35 sponse for some suitable energy. Since the corrected data 

is then logarithmically linear, that data can be utilized 
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by the reconstruction algorithm to produce a spatially- 
consistent image that is independent of patient size. This 
correction may also be implemented by assuming an average 
composition of tissue and bone instead of water for the 
5 various degrees of attenuation. 

Finite Dimensions of Scanning Aperture 
The qualities of resolution and image sharpness are 
closely associated with aperture size. Theoretical devel- 
opment of the reconstruction algorithms is based upon an 

10 infinite amount of infinitely thin transmission data which, 

in practice, is approximated by a finite number of trans- 
mission measurements of an x-ray beam of finite dimensions. 
The minimum aperture size is limited by photon statistics 
for a given x-ray output, geometry and scan speed. 

15 The slit height, perpendicular to the linear motion of 

the scanner, acts to make the resultant reconstructed 
coefficient a type of average coefficient over many thin 
transverse planes. This average is not inconsistent with 
the assumption of infinitely thin rays since the theoreti- 

20 cal development of these algorithms is limited to two dime- 

nsions. Distortions due to this vertical smearing are 
somewhat minimized because of the homogeneity of the human 
body over short vertical distances. In some respects, this 
smearing is advantageous in that a larger volume is consid- 

25 ered in each cross-sectional slice, hence fewer slices need 

be reconstructed to include the entire volume of interest. 
However, in the two-dimensional reconstruction, the resul- 
tant coefficient need not pertain to only one type of 
tissue, especially for small structures and near bound- 

30 aries. This necessity is particularly important in relat- 

ing these coefficients to effective atomic number, density, 
chemical composition, or specific tumor types. 

The slit width, parallel to the linear motion of the 
scanner, introduces a compromise between the algorithm 

35 development and practical measurement. The approximation 

that the average intensity transmitted along the width of 
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the aperture is the same as the relative intensity along a 
central ray is partially responsible for the lack of 
sharpness noted along many boundaries and can significantly 
degrade resolution. If the slit width is less than the 
5 linear increment between samples, computer enhancement of 

the data is not feasible. However, if the slit width 
exceeds this increment, then information is available from 
which the intensity for, an aperture of a width approaching 
that of the linear increment can be calculated in accor- 

XO dance with the present invention. This technique is 

comparable to reconvolving the point spread function from 
other imaging devices, except the correction is applied to 
the measured data before an image is formed rather than as 
a modulation transfer function enhancement done as a post 

15 processing procedure on the final reconstructed image. The 

deconvolution is along discrete steps, the measured data, 
and not a continuum, thereforeading to the limit of data 
resolution, the linear increment. 

Time Response of Detector System 

20 The time response is an important consideration in the 

choice of detector systems. If signal decay due to an 
impulse of radiation on a detector is slow, then it would 
appear from the measurements alone that some radiation was 
still incident on the detector some short time after 

25 radiation exposure. Likewise, a particular measurement 

during a scan may be partially due to radiation incident on 
the detector from some time prior to that measurement. In 
the case of the linear scanner considered here, this, 
temporal smearing due to the detector response may be 

30 related to a type of spatial response of the detector by 

the time spacing of the measurements or scan speed. In 
this way, the time response of the detector is considered 
similar to or as part of the aperture transmission function 
and is corrected for in the same manner. 

35 Applicability to Fan Beam Scanners 
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For fan beam scanners, many of these correction tech- 
niques still apply. Furthermore, the temporal response of 
the detector is related to measurements made with the same 
detector as a function of angle (time) . The aperture 
itself is less likely to overlap than in the case of the 
linear scanner, however, the effective aperture due to 
cross-talk and patient scatter may also be handled by these 
preprocessing techniques. Modifications applicable to 
patient scatter may include making a correction to the 
measurement of interest as a function of intensity attenu- 
ated (including scatter) rather than intensity transmitted 
to a nearby detector along with the distance of this 
detector from the detector of interest. In any case, this 
type of scatter correction is just an approximation. The 
corrections pertaining to geometrical errors, instability 
of the x-ray source and the polychromaticity of the x-ray 
beam are applicable to either type of scanner. 

Some of these techniques are known to be in commercial 
use, especially in regard to polychromatic correction. 
Hardware approaches to some of the problems are also 
common, such as careful construction of the scanner and use 
of patient supports and restraints to minimize geometrical 
error, reference detector methods to compensate for x-ray 
source variation, temperature compensating amplifiers to 
reduce drift, increasing x-ray filtration to reduce poly- 
chromatic effect, increasing x-ray output to allow for 
smaller apertures, and the selection of detectors to reduce 
afterglow. Even so, the present method for reducing theae 
factors which degrade the quality of the data is beneficial 
and in some cases, absolutely necessary. 

In spite of these many improvements and alternative 
approaches, there remains a need for improved methods for 
reducing the effect of these error factors on a back 
projected image, as well as for actually back projecting 
the image of the target object. There is also a need for 
an apparatus for implementing these methods, and especially 
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an apparatus which is capable of being mounted in combina- 
tion with an existing x-ray simulator without changing the 
function of that apparatus. Such an attachment makes 
possible a highly desirable economy and flexibility of 
5 application, especially when combined with the improved 

methods of correcting for the effect of the above-listed 
error factors and back projecting the image of the target 
object. It is, therefqre, an object of the' present inven- 
tion to provide such an apparatus and such a method. The 

10 method utilizes a combination of Hilbert and Fourier 

transforms to produce results which approach the theoreti- 
cal accuracy which the Radon method might produce with 
perfect data and an infinite detector array, and may be 
used for either parallel or fan beam geometries. 

15 Another object of the present invention is to provide 

a method and apparatus which improves the resolution of a 
CT image by providing an improved method for back project- 
ing that image. 

Another object of the present invention is to provide 

20 an improved method for correcting for the errors in the 

calculated back projected image caused by off center or 
eccentric rotation of the source of the x-ray beam and/or 
the detector array. 

Another object of the present invention is to provide 

25 a single x-ray apparatus capable of performing both CT 

scanning and panoramic imaging. 

Other objects, and the advantages, of the present 
invention will be apparent to those skilled in the art from* 
the following description of the presently preferred 

30 embodiment thereof. 

SUMMARY OF THE INVENTION 
These objects are achieved by providing an attachment 
for back projecting a CT scan using an existing x-ray 
35 simulator comprising a linear detector array, the mounting 

of which is adapted for mounting to the film holder of an 
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x-ray simulator normal to the central axis of the fan beam 
produced by the simulator without altering the original 
function of the x-ray source of the simulator, which 
comprises a plurality of individual elements. Each indi- 
5 vidual detector element produces an output signal having an 

amplitude proportional to the energy intensity of an x-ray 
incident thereon, the intensity of an incident x-ray being 
proportional to the density of the target through which the 
incident x-ray passes before striking the individual 

10 detector element. Also provided is means for collecting 

the output signals of each individual detector element at 
a plurality of incremental angles as said detector array 
and the source of the fan beam are rotated around the 
target object and translating the signals into a back 

15 projected computer tomographic scan of the target object. 

Translation is accomplished by correcting for the 
spread of the x-rays in the fan beam and for the differenc- 
es in the intensity of the x-rays comprising the fan beam 
depending on the position of the individual detector 

20 element in the detector array relative to the central axis 

of the fan beam. Translation also involves scaling the 
output signal to account both for the relative distance 
from the source of the fan beam to the detector array and 
for the distance from the fan beam source to the target 

25 object. The attachment also includes means for correcting 

the scaled signal for the magnification resulting from the 
spread of the fan beam and for transforming the output 
signal from each detector element into an output signal- 
representing the output signal that would have been pro- 

30 duced by each detector element had the incident x-ray 

originated from a parallel beam instead of a fan beam. 
Finally, translation involves converting the transformed 
output signal at each incremental angle of the detector 
array into a gray scale value for a picture element having 

35 a specific set of coordinates relative to the coordinates 

of the detector array and outputting the gray scale value 
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to an appropriate display means for displaying the tomo- 
graphic scan. 

Also provided is a method of producing a computed 
tomographic scan of a target object with an x-ray simulator 
comprising the steps of projecting an x-ray beam through a 
target object, detecting the x-ray beam incident on a 
detector positioned on the other side of the target from 
the source of the x-ray beam, defining polar and cartesian 
coordinate systems to describe the geometry of the x-ray 
beam, target object, and detector and locating the position 
of the source of the x-ray beam, the center of rotation 
axis, and the detector on the polar and cartesian coordi- 
nate systems for each incremental angle as the source of 
the x-ray beam and the detector rotate relative to the 
center axis of the target object. The output signal from 
each detector is then transformed into an output signal 
representing the output signal which would have been 
produced by that detector had the incident x-ray originated 
from a parallel beam source rather than a fan beam x-ray 
source and the transformed output signal is converted at 
each incremental angle of the detector into a gray scale 
value for a picture element having a specific position on 

the cartesian coordinate system. The gray scale values 

are then output to an appropriate display device. 

BRIEF DESCRIPTION OF THE DRAWINGS 
Figure 1 is a schematic, sectional view of an x-ray 

simulator having an apparatus constructed in accordance 

with the present invention mounted thereto. 

Figure 2 is a schematic of the component parts of the 

apparatus of the present invention and the processing of 

the data by those parts. 

Figure 3 is a flow chart of the processing of the data 

in accordance with the method of the present invention. 

Figure 4 is a schematic representation of the geometry 

of the coordinate system in which the x-ray source, target, 

and detector array of the present invention are located. 
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Figure 5 is a schematic representation of the geometry 
of the back projected image computed in accordance with the 
method of the present invention. 

Figure 6 is a graphical representation of an exaggera- 
5 tion of a measured center of rotation shift for an x-ray 

simulator rotating 360°. The off-center shift Sj is plotted 
as a radius from the center versus the angle in degrees. 

Figure 7A is a graphical representation of a method 
for measuring the center of rotation shift as a function of 
10 an angle j ; Figure 7B is a schematic representation of the 

type of error in data acquisition and back projection, 
which results in distortion and blurring of the back 
projected image of a phantom, as a result of the superposi- 
tion of translational and rotational movement of the x-ray 
15 source of a panoramic dental x-ray device with respect to 

the phantom. For purposes of clarity, only two rays P x 
(shifted to P 1 + 5) and P 2 (shifted to P 2 + 5) are shown as 
having been shifted by this eccentric rotation. 

Figure 8 is a schematic representation of the diver- 
20 gence of a ray emanating from a point source and the 

different magnification factors that result as the ray 
passes through the object. 

Figure 9 is a schematic representation depicting the 
circle of image reconstruction showing the grid defining 
25 picture elements (pixels) and the area enclosed by two 

superimposed projections of the fan beam. 

Figures 10A and 10B show alternative methods using the 
spacing element Sr which is utilized during back projection 
as either a constant (Fig. 10B) or variable value (Fig. 
30 10A) . 

Figure 11 is a representation of the rectangular pixel 
grid containing the circle of image reconstruction and 
showing the sum of the back projected geometrical factors 
at each pixel in the 16 x 16 grid. 
35 Figures 12A and 12B illustrate the filters used for 

processing the output signals of the detector elements and 
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the effect of each type of filter on the convolved projec- 
tion data, respectively. 

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 
Referring to Fig. 1, there is shown a schematic repre- 
sentation of an x-ray simulator, indicated generally at 
reference numeral 10. X-ray simulator 10 is generally com- 
prised of the x-ray tube 12 used for radiation therapy 
simulation which produces a fan beam, indicated at refer- 
ence numeral 14, incident upon a target object 16 posi- 
tioned on table 13. X-ray beam 14 is also incident upon 
the detector array (not shown) of the apparatus 20 of the 
present invention. The apparatus 20 is comprised of the 
linear detector array 26 (see Fig, 2) , the mounting of 
which is adapted for fitting into, mounting next to, or 
covering the film holder 22 of x-ray simulator 10 normal to 
the central axis of the fan beam 14 produced by the source 
12 without altering the original function of simulator 10. 
The film holder 22, or cradle, is mounted on top of the 
image intensifier tube 24. The x-ray beam is collimated by 
using the shutters or collimators (not shown) provided by 
the x-ray therapy simulator 10. The beam 14 is preferably 
shaped with a width of approximately one centimeter and a 
spread determined by the length of the detector array 26 of 
the apparatus 20 of the present invention. If the project- 
ed shadow of the target object 16 to be scanned is less 
than the maximum width of the detector array 26, the spread 
of the beam 14 is decreased so as to just encompass the 
diameter of the target object 16. 

Fig. 1 is shown as a schematic representation of a 
conventional x-ray simulator which rotates around a target 
object (as will be explained below) in a 360° arc for 
reconstructing an image. However, those skilled in the art 
who have the benefit of this disclosure will recognize that 
reconstruction can also be accomplished using projections 
from an arc of rotation that is less than 360°; for exam- 
ple, by using a dental panoramic x-ray unit which rotates 
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through an arc of approximately 230-240°. It is specifi- 
cally contemplated that the apparatus 2 0 can be mounted to, 
and the method of the present invention can be accomplished 
in connection with such dental panoramic units as, for 
5 instance, the AUTOPAN (Belmont Equipment Co. , Summerset, 

N.J.)/ GX-PAN or PANELIPSE II (Gendex Corp., Milwaukee, 
WI) , PANOREX II (Keystone Corp., Neptune, N.J.), or, in a 
presently preferred embodiment, the Siemens Medical Systems 
(Charlotte, N.C.) OP 10 dental x-ray unit. 

10 Referring now to Fig. 2, there is shown a schematic 

representation of an exemplary apparatus used in connection 
with the method of the present invention. The detector 
array 26 utilizes an image intensifying screen (not shown 
in Fig. 2) placed over a linear photodiode array of detec- 

15 tors. In a presently preferred embodiment, the photodiode 

array is comprised of eight to sixteen (or more) light 
sensitive modules (not shown) , each composed of a silicone 
chip with 127 photodiodes (not shown) such as the Model 
THX-1086 (Thomson Electron Tubes Division, CSF Thomson et 

20 Cie (Paris, France) . The photodiodes are used in a capaci- 

tive storage mode (with reverse bias) so that the light 
emitted by the intensifying screen discharges the photodi- 
ode capacitors. The scintillator material used is ytterbi- 
um-gadolinium oxide, a screen material commonly used in 

25 diagnostic radiology. Other scintillator materials may 

likewise be used to advantage, and the screen may be made 
of different thickness to improve photon collection effi- 
ciency. 

An external sync trigger circuit 28 is provided for 
30 interfacing the x-ray simulator 10 with an external sync 

circuit 3 0 to trigger the timing circuitry of the apparatus 
20 of the present invention when the x-ray tube 12 starts 
to emit the fan beam 14. External sync circuit 30 provides 
the synchronization signal to the preprocessor 32 and 
35 determines the integration period of the detector array 26, 

The external sync circuit 30 consists of an input (not 
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shown) to provide software control of offset and gain 
calibration reset signals to the preprocessor 32 , a contin- 
uous clock circuit (not shown) to provide the synchro- 
nization and integration period for the detector array 26 
during gain calibrations, and a variable clock circuit 
(also not shown) which utilizes an EPROM to provide the 
synchronization and integration period for detector array 
26 during data acquisition. Preprocessor 32 receives the 
analog signal from the detector array 26 and, under control 
of the external sync circuitry 30, digitizes and outputs 
the data to the microcomputer interface circuit 34. 

Computer interface circuit 34 converts the data 
received from the preprocessor 32 to a format compatible 
with the input-output (I/O) board 36 of the particular 
microcomputer utilized for processing the data as described 
below. The microcomputer I/O board 36, after receiving the 
data from microcomputer interface circuit 34, inputs the 
data to the memory of the microprocessor 38, which is 
preferably a Motorola 6803 0 or Intel 80386-based computer. 
In a presently preferred embodiment, an Apple Mac II Plus 
(tm) with eight megabyte memory, large hard disk, 80-140 
megabytes, floppy disk and CRT displaying a 480 x 512 pixel 
image (or, alternatively, a special CRT to display a 512 x 
512 image) having a 25 and preferably 3 0 megahertz internal 
clock is used to advantage. In addition, the output of the 
microcomputer I/O board 36 provides the software control 
signals to the external sync circuit 30 for offset and gain 
calibrations, the continuous clock and variable clock 
circuits. The microcomputer 38 performs the data manipula- 
tion routines described below to produce the CT image, 
outputting the processed data to an appropriate high 
resolution monitor 40 which provides a gray scale or color 
display of the reconstructed CT image, the film recorder 
(matrix camera) 42 which provides a hard copy output (film) 
of the image on the high resolution monitor 40, and/or the 
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large capacity storage device 44 which stores the image 
data as well as patient information . 

A general schematic outline of the data manipulation 
process is set out in Fig. 3. By the term "data" it will 
be understood that reference is made to the processed 
output signals of each individual detector or photodiode 
comprising detector array 26, which is proportional to the 
intensity of the x-ray ^incident thereon. Preprocessor 32 
integrates the output signal of each photodiode and reads 
out the resultant signal upon external command controlled 
by microcomputer 38. In so doing, the data is corrected, 
smoothed, and filtered as shown in the step represented by 
box 46 in Fig. 3. This filtering step is required because 
the individual detector elements, or photodiodes', become 
defective and may malfunction either partially or complete- 
ly. A correction algorithm is provided which identifies 
those faulty detector elements through comparison with the 
neighboring detector elements (the narrow separation 
between detector elements utilized allows the generaliza- 
tion that the response for adjacent detector elements 
should be relatively close to each other) , and replaces the 
bad data with a corrected, interpolated value. A smoothing 
algorithm is provided for extending this method by setting 
a limit, or preset selection criterion, on the variation of 
adjacent individual detector elements and interpolating 
response values for those detector elements whose variation 
from the neighboring detector elements exceeds that prese- 
lected criterion. Recognition of when an individual, 
detector element malfunctions is determined experimentally 
by exposure to x-rays both in air and through a phantom 
target (not shown) and the history of the detector array. 

Characterizing each step in more detail, the data is 
first smoothed in accordance with the formula 

P(2) = (P(l) + P(3))/2 [1] 
where P(l), P(2) , and P(3) represent the signal outputs 
from each of three successive detector elements in detector 



WO 93/00649 



PCT/US92/05276 



array 26. The smoothing is extended by selecting a maximum 
variation between adjacent detectors based upon the detec- 
tor response history. For example, a cutoff limit of ± 0.2 
or 0.3 times the adjacent value may be selected so that if 
the reading from a particular detector element varies from 
its neighbors by more than 0.2 or 0.3 times the output of 
the neighboring elements, that output signal is replaced by 
the average of the output signal from the neighboring 
elements on either side of the bad element. In the event 
that several detectors in a row have malfunctioned, a 
scaled averaging method is used to correct the data. The 
difference between the last good output signal at one end 
of the malfunctioning row of detector elements and the 
first good output signal at the other end is divided by the 
number of detector elements malfunctioning and the resul- 
tant guotient is used as a constant which is sequentially 
and iteratively added to the last good output signal to 
provide a replacement for the next output. The correction 
is added to this value for the next reading and so on, up 
to the first good output signal. For example, given the 
following output signals, 110, 125, 255, 255, 255, 255, 
140, 141, the correction is applied as follows: 

C = (140 - 125)/5 = 3 [2] 
In this example, there are five intervals from the last 
good output signal at one end of the malfunctioning detec- 
tor elements and the first good output signal on the other 
end of the malfunctioning detector elements. The corrected 
output data is 110, 125, 128, 131, 134, 137, 140, 141. 

Other corrections account for variations in x-ray 
output. These variations can be caused by current and/or 
voltage variations. Metering circuits on the x-ray simula- 
tor 10 are used to monitor these parameters, and the output 
of those circuits is used as an input to processor 32. 
Current variations are corrected in linear fashion: if, for 
instance, the current drops from 5 mA to 4 mA, the detector 
elements are corrected as (5/4) x P(i) . Voltage correc- 
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tions are non-linear. The correction factor is formed as 
a power function of the voltage, i.e., V N , where N is a 
number in the range 1.5 < N < 3. Nis determined experi- 
mentally and is idiosyncratic to the system being con- 
5 trolled. Suppose the voltage drops from 120V to 110V and, 

for that range, it has been determined that N = 2. The 
correction is made as follows: 

P(i) corr - [120/110] 2 x P(i) [3] 
Additional corrections are performed using interpola- 

10 tion procedures to optimize accuracy in the calculated 

intensities used in back projection when correlating the 
detector in the detector array -26 that intercepts a back 
projected, or parallel beam, ray 50 (Fig. 4) which passes 
closest to the point in the target object 16 being recon- 

15 structed. The coordinate system in the back projection 

process is defined from i = 0 at the first (-x) detector to 
N D - 1, where N D is the number of detectors in the array 26. 
For example, using a 2048 detector array 26 with a detector 
spacing of x = 0.045 cm, the total length D L of the line of 

20 detectors in detector array 26 is 

D L = (2048) x 0.045 cm = 92.16 cm. [4] 
The factor 2 048 is used because the detector line is based 
on the spacing of the individual detectors and there are 
2047 spaces between the 2048 individual detectors in a 

25 presently preferred embodiment of detector array 26. The 

center of the detector line is taken as zero. The x 
coordinate along this line then varies from -92.16/2 to 
92.16/2, and the x coordinate of the i th detector is 

x L = -46.08 cm + (i + 0.5) x 0.045 cm, [5] 

30 e.g., for the 1135 th detector, 

X 1135 = -46.08 + (1135 + 0.5) x 0.045 = 5.0175 cm. 
The coordinate on the detector line where the ray 
passes closest to the reconstruction target point is desig- 
nated L bp (back projected) . Once this coordinate is 

35 determined, it is possible to calculate which individual 

detector will receive the ray passing through L bp . Let 
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N bp = Cf^p + 46.08)/<fx) + 0.5] [6] 
and N bp = INT[N bp ], 

where "INT[]" indicates the integer of the number in 
brackets and the factor 0.5 is an end correction such that 
5 in the limit where is approximately equal to -46.08, the 

formula will return the integer value of one (1) . 

Round off error in calculating this integer is mini- 
mized by linear interpolation. The fractional part of N bp 
is 

10 N bp - [N bp ] = f bp . [7] 

In the simplest round off, i.e., if f bp < 0.5, N bp = [N bp ] . 
If f bp > 0.5, N bp = [N bp ] +1. For example, if N bp = 756.528, 
the detector output taken would be PI (757) . To be more 
exact, interpolation yields 

15 PI = (1 - f bp ) Pl([N bp ]) + f bp Pl([N bp ] + 1). . 

In the numerical example, PI = (1 - .528) PI (756) + .528 
Pl(757) . 

The term "smoothing" (steps 46 and/or 72 of Fig. 3) 
refers to the filtering of the data by both smoothing and 
20 ripple filtering. Smoothing is accomplished using a set of 

reconstructed values which may be represented as 

M(l,l) M(l,2) 11(1,3) M(1,N) 

M(2,l) M(2,2) 11(2,3) M(2,N) 

M(3,l) M(3,2) M(3,3) M(3,N) 

25 

M(N,1) M(N,2) M(N,3) M(N,N) [8]. 

Here, N = 2 k where k = an integer, i.e., N = 128, 256, 512, 

30 etc. For four point smoothing for nearest neighbor (j 

denotes the row and i denotes the column) : 
M(j,i) = [M(j-l,i)+M(j+l,i)+M(j,i-l)+M(j,i+l)]/4 [9] 
The top, bottom, and adjacent left and right neighbors are 
used with equal weight. For five point smoothing for 

35 nearest neighbor, 

M(j,i) = [M(j,i)+M(j-l,i)+M(j+l,i)+M(j,i-l)+M(j,i+l)]/5 [10] 
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M(j,i) is added in along with the four nearest neighbors. 
For weighted neighbor smoothing with the central value 
given the greatest weight, nine point smoothing is used: 
M'(j,i) = [lxM(j-l, i-1) + 2xM(j-l,i) + lxM(j-l,i+l) 

+ 2xM(j,i-l) + 4xM(j,i) + 2xM(j,i+l) [11] 
+ lxM(j+l,i-l) + 2xM(j+l,i) + lxM(j+l,i+l) ]/16. 
This weighting system can be varied with different weights 
applied to different locations depending upon the circum- 
stances as known to those skilled in the art. 

Ripple filters are used to eliminate the wave intro- 
duced into the reconstruction by employing Fourier- methods. 
The first step in application of such a filter is to obtain 
x-ray transmission values for a water (uniform target 
density) phantom. With known reconstructed CT values, 
these true values are compared to observed CT values. The 
ripple filter factor for each picture element is given by 
the true value /observed value. To achieve correction, the 
value for each picture element is multiplied by the filter 
factor. There are many other different types of filters 
which may be used to advantage depending upon the circum- 
stances as known to those skilled in the art who have 
benefit of this disclosure. Such filters might include, 
for instance, edge enhancement filters, high frequency 
suppression, bone enhancement filters, soft tissue enhance- 
ment filters, and so on. 

The smoothed, corrected data from step 46 is convolved 
at step 47 and then back projected at step 48. Convolution 
38 and back projection 48 are accomplished in accordance 
with an algorithm which is derived as follows. Referring 
to Fig. 4, there is shown the geometry for the fan beam 14, 
the parallel beam 50 and the coordinate system used to 
describe the algorithm of the present invention. In this 
system, the origin is established at the nominal center of 
rotation 56 with the knowledge that the nominal center may 
vary during the rotation. The x-axis 52 is on a line 54 
through the center of rotation 56 so that the +x axis will 
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coincide with the left side of the target object 16 and 
perpendicular to the y-axis 58 which passes through the 
center of rotation 56 on a line from the x-ray source 12 to 
the center of the detector array 26 (which is parallel to 
the x-axis 52) , with the +y axis oriented toward the front 
of target object 16, As shown in stick figure in Fig, 4, 
if target object 16 is a patient, the +z axis is* oriented 
toward the patient's {lead; the coordinate system, when 
defined in this manner, is a right-handed system. The x- 
ray source 12 , is assumed to rotate counterclockwise about 
the z axis with the angle of the x-ray defined from the +x 
axis, and the starting point of rotation is assumed to be 
90° , or the x-ray target coinciding with the +y axis. In 
this fashion, images will be reconstructed that show the 
right side of the suject on the left side of the image and 
reflect the viewpoint of looking at the image from the foot 
of the patient toward the head. 

The distance from the x-ray source 12 to the center of 
rotation 56 is D, reference numeral 60. The separation 
between the center of rotation 56 and the plane of the 
detector array 2 6 is a distance V, reference numeral 62. 
The x-ray source 12 and detector array 20 are incrementally 
rotated about the center of rotation 56 through angles T as 
shown at reference numeral 64 in Fig. 5. The increments 
may be equal in size or unequal. When equal increments are 
used, T = 180°/N or T = 360°/N, where N is the number of 
increments, and the system may be rotated either through a 
semi-circle or the entire arc of a circle. Rotation, 
through the arc defining a semicircle is commonly used in 
dental panoramic x-rays, and the coordinate system defined 
above is equally applicable in that instance. In such 
apparatus, the center of rotation 56 is not fixed such that 
both D and V vary at each angular increment. As set out 
below, the change in the location of the center of rotation 
56 is measured and, using the method of the present inven- 
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tion, the back projection is corrected for the resulting 
known error. 

For equal increments, the index, j, will nominally run 
from 0 to 3 59 or from 0 to 719 depending upon whether the 
5 projections are taken at one degree or 0,5 degree inter- 

vals. The number of angles used for acquisition of data is 
limited by the ability of the data collection system and 
the computer system to acquire and process the given volume 
of input data. The increment between angles Tj and Tj +1 is 
10 denoted by ST and is given by 360/N or 720/N. 

The projection angle, T j , indicated at reference 
numeral 66, is given, in the equal increment case, by 



Tj = 90° + j£T where 0 < j < N - 1. [12] 
The x-ray source coordinates for each j are given by 

15 X g (j) = D g (j)cos(Tj) . [13] 

y s (j) = D s (j)sin(Tj) . [14] 



Since the origin of coordinates is at the center 56 of the 
target 16, the nominal distance of the x-ray source to the 
origin is 

20 D g (j) = [x s 2 (j) + y s 2 (j)] 1/2 

The coordinates, x g (j) and y s (j) must be known for all 
angles Tj . The x-ray source and detector array are assumed 
to be rigidly connected and co-moving in the rotation about 
the target 16. As described above, the angle T 0 = 90° 

25 defines the initial configuration of the system before 

rotation, e.g., at j = 0. At this angle, T 0 , the coordi- 
nates of each of the detectors in detector array 26 are: 

x d (0,i) = - D L /2 + (i + 0.5) Si [15j 
y d (0,i) = - V [16] 

30 where Si = D L/ N d' d l ^ s the l en< ? th of the detector line, and 

N d is the number of detectors in array 26 where 0 < i < N d . 
Equations [15] and [16] describe the coordinates for non- 
eccentric motion; as described above, when motion is 
eccentric, both x and y can be measured and input or solved 

35 for iteratively. The 0th detector is considered the left 

most detector on the -y axis, and the spacing is actually 
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the distance between the centers of each detector in array 
26. The factor 0.5 in equation [15] serves to define the 
end point of the ray in the center of the detector element. 
For Tj > 90° , the coordinates of each detector in detector 
5 array 26 are given by: 

X d (j,i) = x d (0,i) cos (j*T) - y d (0,i) sin (j<ST) [17] 
Y d (J,±) - * d (0,i) sin (j*T) + y d (0 ri ) cos (j*T) . [18] 
In Figure 4, the parallel rays 50 are constructed by back 
projection from the detector array 26. The intersections 
10 of the parallel 50 and fan beam 14 rays define the circle 

of reconstruction 68 . 

In defining the circle of reconstruction 68, the next 
step 70 (see Fig. 3) in manipulation of the data is accom- 
plished by microprocessor 38 using the following steps. 
15 The diameter D rc of the reconstruction circle 68 along the 

line 54 is given by the scaling relationship 

D rc /D = D L /W. [19] 
The radius of reconstruction circle 68 is R rc = D rc /2 and W 
= D + V defines the distance from the x-ray source 12 to 
20 the center of the detector array 26. The Jacobian of the 

scaling transformation from the detector array 26 to the 
line 54 along the diameter of reconstruction circle 68 

M x = D/W [20] 
The scaling factor M 2 allows transformation of the data 
25 taken along the detector line to appropriate values along 

the parallel line 54 through the center of rotation 56 and 
lying along a diameter of the reconstruction circle 68. 

The distance between successive back projected par- ♦ 
allel beams or detector rays 50 at the level of the diame- 
30 ter of the reconstruction circle 68 is the same as the 

individual detector aperture projected back to the center 
of rotation as a length 

Si' = Si M x [21] 
The scaling factor may be written explicitly for each row 
35 ioo (see Fig. 8) of pixels in reconstruction circle 68. 

An index c is defined such that c = 0 indicates the top row 
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100 of picture elements in the reconstruction circle 68 and 
c = N p - 1, e.g., c 15 at reference numeral 104 in Fig. 8, 
indicates the bottom row. Similarly, the index b used in 
equation [38] below indicates the columns 102 of pixels in 
5 the reconstruction circle 68. Here, b = 0 indicates the 

left most column. The magnification M x defined for the row 
through the center of the reconstruction circle 68 is 
generalized as * 

Magp(c) - [D - y(c)]/W [22] 

10 where y(c) - D - R rc + (c + 0.5) D py . 

D py is the height 106 of the picture element of row c and 
D px is the width 108 of the picture element of row c. The 
row of pixels lying on the diameter of the reconstruction 
circle 68 therefore have width and height given by 

15 D px or D py = D rc /N p [23] 

where N p 2 = number of reconstruction pixels; i.e., 256, 512, 
etc. 

To normalize the projected data from the incident 
beam, fan beam projection data is taken for x-ray transmis- 
20 sion through air and through target objects. Let 

p(j,i) = log (r a(j/i) /r p(jfi) ) [24] 
where r a( j #i) is the ith detector reading at angle Tj in air 
and r p( j /i} is the ith detector reading at angle Tj for x- 
rays transmitted through the patient or other target 16. 
25 The ranges of the integer indices are as defined above. 

Projection data for each angle Tj and detector i has 
been defined as p(j,i). This function is re-written in 
terms of the spatial variable s and angle <f> in the reconn 
struction circle 68 as p(s,tf>) . The one dimensional Fourier 
30 transform of p(s,0) at a given angle <t> is written as 

irp]{S,Q) = J" p(s,4>) exp[-27riSs] ds [25 



In the discrete form, for each angle Tj, the one dimensional 
Fourier transform of the projection data is 
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N n -1 

[5" p] (S,*) = I p(k,0) exp[-27riss(k) ] [26] 
k=o 

where s(k) = -R rc + (k + 0.5) 6s. Recall Ss = 2R rc /Njj and Ss 
- the increment between adjacent rays in the spatial 
domain, where s is the spatial domain variable, S is the 
frequency domain variable, i = V-l in this equation, and s 
= 1/S. In this application, s is expressed along the line 
R (65 in Fig. 5) in the reconstruction circle 68, and R is 
displaced on angle <p , measured counterclockwise from the +x 
axis. Integration is taken over the entire real domain 
where -« < S < 

The one dimensional Fourier transform of the function 
p(s,p), written only as a function of s for a given 0, is 
represented in the continuous notation as 

[&~p]S = j m pis) exp[-2TtiSs]ds [27 

It is necessary to perform a two dimensional inverse 
Fourier transform, to obtain the linear absorption coeffi- 
cients, ji(x,y), leading to the calculation of target 
optical densities. Coordinates (x,y) are used in the 
Cartesian spatial domain containing the reconstruction 
circle 68 but defined by the axis along the center of the 
x-ray beam and the orthogonal detector array. The equiva- 
lent coordinates in the frequency domain are (X,Y) . The 
projection at angle Tj for detector element i, p(j,i), is 
expressed as a function of the new coordinates. Thus, 
p(j,i) is written as p(x,y) in the spatial domain and* 
P(X,Y) in the frequency domain. In practice, the one 
dimensional transformation along the line s is equivalent 
to the two-dimensional transformation for points on the 
same line at the angle <p. Thus, the inverse Fourier 
transform, ? " x in this notation, is 



ft 
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where X +M , Y+ M are the upper and lower bounds where P(X,Y) 
exists . 



In the polar coordinates (S,*), this equation becomes 
[7 - 1 P](s,0) = 



P(S, <&) exp [2xi(xS cos$+ yS sin$) 3 \s\dSd$ [29] 



Note the appearance of |S|, the ramp function, as a result 
of applying the Jacobian to transform the coordinate 
systems. The ramp function acts as a filter* Here, the 
frequency is determined by the sampling theorem, which 
states that where a function h(s) is defined over the range 
2R rc (diameter of the reconstruction circle 68) then the 
Fourier transform of h(s) , H(S), is fully described by 
points <SS = l/2R rc apart (the Nyguist frequency spacing) . 
Conversely, if the range of interest of H(S) is 2S M , then 
h(s) may be sampled at intervals not greater than s = 1/2S M . 
Thus, in equations [29] and [30], S = 0, S lf S 2 . . . S M , and 
S M = 1/2SS (the Nyguist frequency) . 

The modified projection ray P1(S,*) may be directly 
obtained by applying the inverse Fourier transform to the 
result obtained in Equation [26] or by using the convolu- 
tion theorem. The inverse transform is performed in the 
spatial domain. The inverse transform of P(X,Y) (written 
symbolically in the frequency domain as P(J,I)) is just the 
projection itself p(j,0). 

The inverse Fourier transform of the ramp function is 
taken by forming an even Fourier series of the ramp func- 
tion in frequency space and using the coefficients of the 
ascending terms as the inverse Fourier transform in the 




The limits in the integrals are rearranged to 
[J - X P] (s,0)= 

fV P{S, <D) exp [27ti UScos$+ySsin$) ] \s\dSd$ 




[30] 
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spatial domain. The sine term does not appear because the 
ramp filter is an even function in frequency space and 
contributions from the negative and positive parts of the 
frequency spectrum account for the factor, two, multiplying 
the summation. In the present method, the first N D terms 
of the infinite, even, Fourier series are used to approxi- 
mate the convolution- The modified projection ray, P1(S,$) 
is written in terms of the original coordinates (J, I) as 

P1(J,I) - rpCjrii) h(|i - i x |) [31] 

i x =o 

where h ( | k | ) denotes the kth term in the even Fourier 
cosine series; its use has the effect of yielding a circu- 
lar convolution even though the summation is expressed 
monotonically and linearly. Note that h(k) is independent 
of angle Tj . 

It is useful to treat the inverse transform of the 
ramp filter .as a continuous function. This treatment is 
accomplished by using the integral rather than the discrete 
transform, or converting the summation in equation [31] to 
its integral analog/ The function, h(k) , is expressed as 

h(k) = 25s[ {1/2)6s \s\cos(sTzkdsS)dS [32] 

where <Ss = the increment between adjacent rays in the 
spatial domain. The integral in equation [32] is a trun- 
cated version of the transform of the absolute filter 
function with cutoffs at ±S M . Consequently, no frequencies, 
outside the ±S M band will be found in the convolved projec- 
tion data. 

This integral is solved to obtain 

h(0) - l/(4*s) for k = 0; [33A] 
h (k) = -l/7r 2 k 2 $S) for k odd. [33B] 
When the projection data is derived from parallel beam 
geometry, Ss is constant. For fan beam geometry, 5s varies 
from pixel to pixel in the reconstruction circle. 
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Using the explicit form of h(k) , the modified projec- 
tion data is solved for using the inverse Fourier two 
dimensional transform by writing the spatial domain convo- 
lution symbolically as 
5 P1QC0) = h(k) * p(k,0) [34] 

to solve the inside integral of the double integral formed 
by the inverse, two dimensional transform (the * symbolizes 
the convolution operation) . The conceptual development 
leading to equation [34] is used to derive the hybrid 
10 reconstruction algorithm used in the method of the present 

invention. Returning to the discrete representation, the 
Fourier inverse transform of the ramp function, |S|, is 
written as: 

Np-1 

15 . Pl(k,0) = p(O,0)/45s - (l/7r 2 5s) I p(a,0)/(k-a) 2 [35] 

a=l 

where if k = odd, a = 0, 2, 4... and if k = even, a = 1,3, 
5 , so |k-a| = odd and (k - a) * 0. More particular- 
ly, if k = odd, 
20 N D -1 

Pl(k,0)={[p(k,0)/(45s)]-[l/(7T 2 5s)j:p(a,5)/(k-Q:) 2 ]} [36A] 

a=0 ,2,4... 

If k = even, 

N D -1 

25 Pl(k,0) = {[p(k,0)/4£s] - [l/2?r 2 (Ss Ep(a0) / (k-a) 2 ] } . [36B] 

a=l, 3 ,5 ... 

These two equations [36A] and [36B] complete the calcula- 
tion of the inner integral in equation [34]. Completing 
the integration of the outer integral is accomplished by 
30 calculating the integral over 0. This process is tradi- 

tionally called back projection. 

Back projection geometry is shown in Fig. 5. The x,. 
y increments are given by: 

6x, Sy = 2R rc /N p [37] 
35 where N p = number of pixels in one dimension, i.e. 32, 64, 

128, etc., and R rc = radius of reconstruction circle 68. In 
similar fashion to equations [15] and [16], 

x = — R rc + (0.5 + b)5x [38] 
y = +R rc - (0.5 + c)*y [39] 
40 and b,c=0, 1, 2, — N p - 1." 
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The linear attenuation coefficient for each x,y or b,c is: 

N-l 

M(b,c) =E Pl(j,i s (x,y)) *Tl [40] 
i 5 =0 

5 where N = number of angles and ST1 is 7T/N for both 180 and 

360 degrees of rotation. Here, i 5 (x,y) is determined at 
each angle, and j is the number of the parallel ray 50 at 
the projection angle j which lies closest to the point 
(x,y). The absorption .coefficient is calculated from the 
10 following relations. 

The angle <p of the ray of interest is defined by 

$ = tan" 1 (y/x) . [41] 
Because <p is used in conjunction with T jf <f> must be defined 
in the same fashion as Tj, or from 0 to 360 degrees. If x 
15 < o and y >0, then 90° < <p < 180°. If x < 0 and y < 0, 

then 180° < <p < 270°. If x > 0 and y 0, then 270°. < <p < 
360°. Let 

L 8 = W r 5 sin{[T j - 0]/[D - r 5 cos (Tj - <*>)]}, [42] 
where r 5 = (x 2 + y 2 )** C 43 ] 

20 and i 5 (x,y) = INT[(L 8 + D L /2)/(i + 0.5)] [44] 

As before, INT[] is the "integer part of". The factor 0.5 
appears in equation [44] because when L 8 = -D L /2, i 5 (x,y) 
must be zero. With L 8 defined, the number of the detector 
corresponding to L 8 is determined using equations [6] and 

25 [44]. For example, suppose L 8 was calculated as 15.463 and 

the length of detector array 26, as it is in a presently 
preferred embodiment, was 69.12 cm with a detector aperture 
of 0.045 cm, then the number of the detector lying on that 
point is given by: 

30 i 5 = INT[((15.463 + 69 . 12/2) /0 . 045) + 0.5] = 1112. 

This calculation indicates that the ray 50 from detector 
1112 passed nearest the point (x,y) that is being recon- 
structed. 

The quantity L 8 is defined differently for the hybrid 
35 reconstruction algorithm of the present invention than for 

the parallel beam geometry. At each angle, Tj, a line is 
passed from the x-ray target position xTj, yTj through the 
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point x,y and from there to an intersection with the x-ray 
detector line. Thus, L 8 is defined as the distance ± from 
the center of the detector line. For each value of x,y, 
there are N back projections, 
5 In practice, a rectangular grid with N p x N p dimensions 

is reconstructed and all values of absorption coefficients 
Mb c lyi n< ? outside the circle of reconstruction 68 are set 
equal to zero. Using the rectangular grid for reconstruc- 
tion greatly simplifies any smoothing routines when they 
10 are utilized. Further, interpolation may be used when a 

* ray 50 does not pass through a point in the back projec- 
tions. 

Referring now to Fig. 6, which is a graphical repre- 
sentation of an exaggeration of a measured center of 

15 rotation shift for an x-ray simulator rotating 3 60° , the 

off-center shift Sj (shown at 7 6c) is plotted as a radius 
from the center of rotation 56 versus the angle in degrees. 
When the x-ray beam which should fall on the central 
detector in the detector array 26 does not pass through the 

20 center of the circle of reconstruction 68, two errors may 

occur. The first error is a linear displacement of the 
detector array 2 6 parallel to the central axis of the fan 
beam 14. The second error is the magnification or demagni- 
fication of the distance between detector arrays (Si) used 

25 in the hybrid reconstruction algorithm described above. 

The latter error occurs when the center of rotation 56 is 
shifted in a direction perpendicular to the line on which 
the linear detector array 26 is positioned (e.g., the 
detector line) . Further, both errors may occur at the same 

30 time. 

Referring to Fig. 7, a simplified method of measuring 
the shift in the center of rotation is illustrated graphi- 
cally. A round rod (represented schematically at reference 
numeral 74 in Fig. 5) which is relatively opaque to x-rays 
35 is placed at the nominal center of rotation shown by line 

78 in Fig. 7A, and data obtained using the method of the 
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present invention is shown as curve 84. An eccentricity in 
the center of rotation causing a shift in the direction 
parallel to the central axis of the projected fan beam 14 
results in projection data for the rod with center shown in 
5 dotted line 80 with projection data 82. The curves 82 and 

84 shown in Fig. 7A represent the reconstructed projection 
data taken over each angle Tj . From this data, the shift 
in the number of detector widths for each angle is mea- 
sured. This shift is used to correct equation [44] as 
10 follows: 

i 5 (x,Y) = INT[(L 8 + D L /2)/i - Sj + 0.5], [45] 
where Sj is equal to the shift in number of detector widths. 
However, the maximum shift in either direction must be 
limited so that the shifted detector ray 50 at the edge of 

15 detector array 26 still passes through the target object 

16. As shown in Fig. 7B, when the central ray (labelled P 2 ) 
of the detector line does not pass through the center 56 of 
the circle of reconstruction 68, an error is introduced 
which blurs or distorts the image of phantom 78 as set out 

20 above. This error results, for instance, from off-center 

rotation of source 12 and detector array 26 through arc 66 
and arc 27 (referring to the center point of array 26) , 
respectively, or, in the case of a dental panoramic x-ray 
unit, from the superposition of rotational and translation- 

25 al rotation, such that the nominal center of rotation 56 is 

shifted to the point 57. The shift is measured and used to 
correct the backpr o j ection or calculated iteratively by 
subjective evaluation of the reconstructed image resulting* 
from each calculation. For example, if ray P 1 at angle j 

30 is detected as ray P x + 5, it is designated P x + 5 and used 

to re-set the projection value at P + 5 and so on for P 2 (P 2 
+ 5) , etc. In the application described above, the magni- 
fication or demagnif ication of Si is ignored. 

In a presently preferred embodiment, however, the 

35 demagnification actually present at each row of pixels 100 

(see Fig. 8) in the reconstruction circle 68 is accounted 
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for in the reconstruction calculation. The actual varia- 
tion of the Jacobian of the transformation is used in the 
method to apply the appropriate value of demagnif ication, 
Magp(c), to insert the exact value of the back projected 
data in the calculation, 

Magp(c) is the Jacobian of the transformation for a 
coordinate system the origin of which is in the detector 
line with points (X',Y'J to a coordinate system located at 
the pixel row c with points (X,Y). For the Jacobian, 



J[(X,Y)/(X',Y')] = 



6X/6X' SX/6Y' 
8Y/6X' 6Y/SY' 



J = 



Magp(c) 
0 



0 

1 



« Magp(c) 



[46] 
[47] 



Thus, at the pixel line c, 

X'd(i) " x d(i) Magp(c), [48] 

where i is the index for the detector line and 0 < i < N D - 
1, starting at the left side of the detector line when j = 
0. This correction for the magnification resulting from 
the spread of the x-ray fan beam 14. can be better visual- 
ized by referring to Figs. 1, 8, and 10. As best shown in 
Figs. 1 and 8, the spread of x-ray beam 14 is such that the 
image of a target object 16 positioned closer to the 
detector line will be smaller than the image of a target 
positioned farther from the detector line. Conversely, in 
the back projection, demagnif ication must be applied to the 
data detected at the detector line. Consequently, the 
width of a detector element Si, when back projected has a 
dimension Si' = SiMagp(c). A related factor occurs in 
treating the spacing element Sr, reference numeral 110 in 
Fig. 10A, or 112 in Fig, 10B. This spacing element, as 
utilized during back projection, is either constant (Fig. 
10B) or variable (Fig. 10A) , the latter being preferred 
because it results in an orthogonal grid of pixels as shown 
in Fig. 8. 



WO 93/00649 PCT/US92/05276 

X-ray magnification of fan beams follows a harmonic 
progression. The magnification factor, Mag(p)c, for each 
row of pixels in the grid 14 in Fig. 8 is given by the 
ratio of the distance from the source 12 (Fig. 1) to the 
5 level of the particular row of pixels divided by the 

distance from souce 12 to detector array 26 (e.g., D+Von 
Fig. 4) . The numerator of this quotient increases by one 
pixel height (D py ) on. each row of pixels, e.g., in a 
harmonic progression. Successive values of the magnifica- 

10 tion factors Magp(c) are shown down the left side of Fig. 

8. These particular factors result from calculation for a 
grid of 16 x 16 pixels for a reconstruction circle 68 with 
the dimensions summarized above, e.g., D = 100 cm, V = 58 
cm, (reference numeral 65 on Fig. 5) = 21.87 cm, and D py 

15 = 2.73 cm with the source at 90° or j = 0. 

In the back projection process, a value Mp as given in 
the grid shown in Fig. 11 represents the sum of all of the 
magnification factors associated with all the rays passing 
through any one pixel is divided into the sum of all the 

20 convolved projection values at that point to give the 

average of the weighted, convolved value at that point: 

c=N -1 

M_ = 1/{(1/N p ) f[l/Magp(c)]}. [49] 

25 T C=0 

This value is multiplied by 7rMagp(c) to obtain the weight- 
ed, linear attenuation coefficient. The coefficients are 
translated into a gray scale using Hounsfield units by 
establishing the intercept point of -1000 corresponding to 

30 the attenuation of the ray due to air. This point is 

represented by the term b in the equation y = mx + b. The 
slope m of this equation is 1000/m w (E) where Mw(E) is the 
measured linear coefficient of water at x-ray energy E and 
x is the weighted, convolved projection values described in 

35 the preceding paragraph. For example, in pixel (7,9), the 

optical density 79.2 (see Fig. 11) is represented as the 
result of the application of projection through thirty 
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angles into 64 detectors. In the preferred embodiment, 
calculation of attenuation coefficients and optical densi- 
ties is accomplished pixel by pixel in each row in the 
direction of the arrow 120 in Fig. 8 first for the row 100 
5 of pixels and then row by row in the direction of arrow 122 

to more efficiently utilize the microcomputer 38. However, 
once the convolution operations are performed, the calcula- 
tion may follow any pixel to pixel order so long as the 
appropriate values of Magp(c) are applied to each pixel. 

10 It is known in the art that convolution filters use a 

quantity <Ss, which is the spacing between rays, that was 
originally defined for parallel beam geometry. However, in 
the method of the present invention, the variable spacing 
of the rays resulting from fan beam geometry is accounted 

15 for without initially using the quantity Sr as is used in 

parallel beam geometry. Instead, this quantity is account- 
ed for, in conjunction with Magp(c), in scaling each row of 
pixels as set out above. 

The ray spacing appropriate for each pixel is account- 

20 ed for in the convolution using combinations of four 

filters, ramp, sine, Cosine and Hamming (the functions of 
each of these four filters are shown in Fig. 12A) . The 
method for using any of these four filters follows that set 
out above for the ramp filter, that is, the filter function 

25 is convolved with the projection data to produce convolved 

projection profiles at each pixel, and is not repeated 
here. The convolution filters are the Fourier, cosine 
series coefficients and are multiplied by the projection 
data to obtain the convolution profiles (see Fig. 12B) . 

30 Additional description of the use of such filters for 

similar purposes can be found in such references as S. Webb 
(Ed.), The Physics of Medical Imaging, Philadelphia: Adam 
Hilger imprint by IOP Publishing Ltd. (1988), Chap. 4. 
Each of these filters has certain desired effects, and 

35 their use is selected to enhance, for instance, bone-soft 

tissue contrast, soft tissue, and so on as described above. 
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Having described the preferred method in accordance 
with the present method will now be described. Initially, 
eccentricities in the rotation of x-ray tube 12 must be 
determined and the individual detectors comprising detector 
array 26 calibrated. A kit (not shown) comprised of 
several rods in a holder, aluminum filter, and a CT phantom 
are provided for that purpose. Placing the center. of the 
phantom for determining motion eccentricities at the 
nominal center of rotation 56, a scan is performed to 
identify any shift in the center of rotation 56. That scan 
also enables verification of the angle encoder output, 
whether that angle encoder (not shown) is integrated with 
the apparatus hardware (see Fig. 2) or output from those x- 
ray simulators which include such an encoder. 

The aluminum filters are then used to determine the 
effective energies and x-ray output for the nominal operat- 
ing voltage of the x-ray simulator. That data is used for 
correcting individual detector output at different voltages 
as at step 46. With a second phantom comprised of differ- 
ent materials, scans are performed at different voltages. 
The results are used to set up a Hounsfield scale for the 
different energies. 

Once the initial measurements have been accomplished 
to load the correction algorithms into microcomputer 38 and 
25 the detectors of detector array 26 calibrated, the 

patient /target object 16 is positioned and the desired 
current and voltage parameters selected. The computer 
program is then initiated to acquire CT data, and on* 
computer ready signal, x-ray output is initialized- first 
30 and the x-ray rotation started. For j=0 to N, where j is 

the angle in degrees and N is the number of measurements 
(360, 720, etc.), the computer is signaled to acquire data 
at approximately one degree intervals (or any other odd or 
even increments as selected by the operator) . For i=0 to 
35 N D -1, where N D is the number of detectors comprising 

detector array 26 and is operator selectable, PP(j,i) is 
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acquired where PP(j,i) represents the output of the analog 
to digital converter (ADC) in microcomputer I/O board 36 
for the detector i at the angle j . Note that the number of 
increments (angles) multiplied by the number of detectors 
will be the number of individual data points and that, for 
instance, where 0 < PP(j,i) < 4095, 4095 corresponds to a 
12 bit ADC, 8191 corresponds to a 13 bit ADC, 16,383 corre- 
sponds to a 14 bit ADC,* and so on. 

Each of the quantities PP(j,i) is then multiplied by 
a correction number obtained from the calibration procedure 
(step 4 6) to adjust the data to the values that would have 
been read if the individual detector were operating per- 
fectly. This correction also remedies the beam hardening 
artifact as discussed above. A reference detector out is 
also used to correct for variation in the output of the x- 
ray unit as a function of the angle j . 

PP(j,i) is the input to the hybrid reconstruction 
algorithm/convolution routine (steps 47 and 48) , and after 
convolution, Fourier artifact filters are applied as 
described above. For c + 0 to b p - 1 and b + 0 to N p - 1, 
the back projection is done (step 38) and then repeated at 
the next b, next c. 

Reconstructed linear attenuation coefficients are then 
converted to gray scale values using the Hounsfield scale 
derived from scanning the phantoms described above in 
connection with the setup procedures. A water phantom is 
scanned and the reconstructed coefficients are assigned 
1000, with air being 0. Other values, such as for bone^ 
are then automatically given a value in this scale. In 
this fashion, the Hounsfield scale is potentially wider 
than that of prior known CT scales, enabling detection of 
variations in CT values not currently available. By way of 
example, if the scan of a water phantom gives a value of 
0.18 per cm for an average linear attenuation coefficient 
of water, 1000 is assigned as the value for water and 
1000/. 18 - 555.55 is obtained as the scale factor for the 
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unit being used. As is known in the art, 0.36 x 5555.55, 
or 2000, is assigned a gray scale of 0 to represent the 
densest scale for this system, or twice the density of 
water, and 2000/255 = 7.843 is the distribution of the gray 
scale number vs. CT numbers: 

0 CT scale » air = 255 gray scale 
1000 CT scale = water = 127 gray scale 
2000 CT scale - twice density of water - 0 gray scale 
The CT values may be displayed along with gray scale 

numbers is desired. 

Although described in terms of the above presently 
preferred embodiments, it is not intended that the scope of 
the invention be limited thereto. Instead, it is intended 
that changes in the specifics set out above 'which do not 
depart from the spirit of the invention described herein be 
included within the scope of the following claims. 
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WHAT IS CLAIMED IS: 

1. A method of back projecting the image of a CT 
scan comprising the steps of: 

(a) exposing a target object to an x-ray fan beam; 

(b) producing a signal representing the intensity of 
a beam of the x-ray incident upon a detector 
after passing through the target object; 

(c) scaling the signal in accordance with the signal 
which would have been produced by the detector 
producing the signal had the x-ray incident 
thereon not passed through the target object; 

(d) correcting the scaled signal for the magnifica- 
tion resulting from the spread of the x-ray fan 
beam; 

(e) convolving the scaled signal to produce a value 
representing the signal which would have been 
produced by the detector had the x-ray fan beam, 
the central axis of which is contiguous and co- 
linear with the central axis of the fan beam; and 

(f) repeating steps (a) - (e) at a plurality of 
incremental angles around the target object to 
produce a back projected image of the target 
object. 

2. The method in accordance with claim 1 of rotating 
the beam source through a circular arc about a center so 
that a plurality of individual exposures are made of the 
target object at angular increments along the arc of 
rotation. 

3. The method in accordance with claim 1 wherein a 
convolved profile is produced for each angular increment by 
a series of steps wherein a discrete Fourier transform is 
performed on each scaled signal and the resulting trans- 
forms are multiplied by increasing values selected from a 
ramp frequency function to provide scaled transforms for 
each angular increment and the inverse Fourier transform of 
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each scaled transform is taken to obtain the value repre- 
senting the parallel beam. 

4. The method in accordance wight claim 3 of using 
the back projected parallel beam established for each 
angular increment to establish a coordinate grid in the 
tomographic plane defined by the intersection of the fan 
beam and the target object, the parallel and fan beams 
intersecting in a locus* of points which defines the diame- 
ter of a reconstruction circle which is coplanar with the 
fan and reconstructed parallel x-ray beams, the coordinate 
grid defined by the spacing between the detectors in a 
linear detector array, and a scaling transformation which 
is the Jacobian of the transformation from the linear 
detector array to each of a plurality of planes through the 
reconstruction circle. 

5. The method in accordance with claim 4 wherein the 
tomographic plane is divided into picture elements, the 
size and shape of which are determined by the coordinate 
grid. 

6. The method in accordance with claim 5 of using 
each back projected parallel beam containing geometrically 
scaled transforms in the back projected geometry to calcu- 
late the linear attenuation coefficient for each picture 
element of the tomographic plane. 

7. The method in accordance with claim 6 of trans- 
lating each said linear attenuation coefficient into an 
assigned shade of gray. 

8. The method in accordance with claim 7 of defining* 
a composite image of the tomographic plane by combining the 
gray shades for all the picture elements. 

9. The method in accordance with claim 1 additional- 
ly comprising measuring the shift in the center of rotation 
for each of the angular increments to correct the values 
for errors caused when the central axis of the fan beam 
does not pass through the center of reconstruction by using 
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a reference rod which is relatively opaque to x-rays place 
at the nominal center of rotation. 

10. The method in accordance with claim 9 wherein the 
rod is used as a target object for yielding projection data 
from which the shift of said center of rotation is derived 
for each angular increment in terms of the number of 
detector widths of each shift. 

11. The method in* accordance with claim 1 additional- 
ly comprising correcting errors of linear displacement of 
the parallel beam relative to the detectors in a linear 
detector array so that a shifted, back projected parallel 
beam at the edge of the detector array will pass through 
the target object. 

12. An attachment for computing a CT image using an 
existing x-ray simulator comprising: 

a linear detector array adapted for mounting to the 
film holder of an x-ray simulator normal to the 
central axis of the fan beam produced by the 
simulator without altering the original function 
of the simulator; 

said detector array comprising a plurality of individ- 
ual detector elements, each individual detector 
element producing an output signal having an 
amplitude proportional to the energy intensity of 
an x-ray incident thereon, the intensity of an 
incident x-ray being proportional to the density 
of the target object through which the incident 
x-ray passes before striking the individual 
detector element; and 

means for translating the output signal of each indi- 
vidual detector element at a plurality of incre- 
mental angles as said detector array and the 
source of the fan beam are rotated around the 
target object into a back projected computer 
tomographic scan of the target object by 
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(1) correcting for the spread of the x-rays in the 
fan beam and the differences in the intensity of 
the x-rays comprising the fan beam depending upon 
the position of the individual detector element 
in the detector array relative to the central 
axis of the fan beam, 

(2) scaling the output signal to account for the 
relative distance from the source of the fan beam 
to the target object and for the distance from 
the source of the fan beam to the detector array, 

(3) correcting the scaled signal for the magnifica- 
tion resulting from the spread of the x-ray fan 
beam, 

(4) transforming the output signal from each 'detector 
element into a signal representing the signal 
that would have been produced by each detector 
element had the incident x-ray been produced by 
a parallel beam instead of a fan beam, 

(5) converting the transformed signal at each incre- 
mental angle into a gray scale value for a pic- 
ture element having a specific set of coordinates 
relative to the coordinates of the detector 
array , and 

(6) outputting the gray scale value to an appropriate 
display means. 
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